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3D topological insulator/ s-wave superconductor heterostructures have been predicted as candi- 
date systems for the observation of Majorana fermions in the presence of superconducting vortices. 
In these systems, Majorana fermions are expected to form at the interface between the topological 
insulator and the superconductor while the bulk plays no role. Yet the bulk of a 3D topological in- 
sulator penetrated by a magnetic flux is not inert and can gap the surface vortex modes destroying 
their Majorana nature. In this work, we demonstrate the circumstances under which only the surface 
physics is important and when the bulk physics plays an important role in the location and energy of 
the Majorana modes. 

PACS numbers: 74.45.+c,74.25.Ha 



I. INTRODUCTION 

Topological quantum computation is one of the most 
active areas of research in condensed matter physics. It 
promises to provide the advantages of quantum com- 
putation such as vast parallelism but with an inherent 
immunity from decoherence. This allows for the for- 
mation of qubits without the need for error correcting 
algorithms 1 . The existence and stability of non-Abelian 
anyons forms the backbone of any architecture for topo- 
logical quantum computation 2 ^. The simplest of these 
excitations is the Majorana fermion. Many diverse 
systems are predicted to harbor these heretofore elu- 
sive excitations including p-wave superconductors 2 ^, 
the v = J fractional quantum Hall state 5 , and cold-atom 
systems^. 

Recently, the search for Majorana fermions has ex- 
panded into the family of materials commonly referred 
to as topological insulators. Generally speaking, topo- 
logical insulators are a class of materials with an insu- 
lating time-reversal invariant bandstructure for which 
strong spin orbit interactions lead to an inversion of the 
band gap at an odd number of time reversed points in 
the Brillouin zone. Topological insulators are differen- 
tiated from other ordinary band insulators by the pres- 
ence of surface states containing Fermi arcs which en- 
capsulate an odd number of Dirac points and are asso- 
ciated with a Berry's phase of n. Normally, such degen- 
eracy points in the bandstructure are easily removed by 
any perturbations, but in the case of topological insula- 
tors the band crossing at the boundaries is protected be- 
cause Kramer's theorem prevents time-reversal invari- 
ant perturbations from opening up a gap in the energy 
spectruirP. In the inceptive work of Fu and Kane 9 , they 
show that coupling s-wave superconduct ors t o 3D time- 
reversal invariant topological insulator d 8 ! 10 ^ via the 
proximity effect may be a potential platform to realize 
these non-Abelian anyons. In particular, Fu and Kane 
show that the surface of a 3D topological insulator - s- 



wave superconductor heterostructure, exhibits many of 
the same properties as a chiral p-wave superconductor 4 
in that the cores of the vortex excitations may harbor 
Majorana fermions. 

Nevertheless, while the analysis presented in Ref . 
considers the gapless surface-state Hamiltonian proxim- 
ity coupled to a superconductor, it ignores the properties 
of the bulk topological insulator. If the bulk were simply 
a trivial insulator it would be inert and no further con- 
siderations would be required. However, it is known 
that topological insulators react to the presence of thin 
flux tubes 13 14 , which can generate a "worm-hole" effect 
that traps low-energy states on the flux tube. This is 
of particular concern as the simplest approach to create 
vortices in a 3D topological insulator - s-wave supercon- 
ductor would be to coat the surface of the 3D topologi- 
cal insulator with a type-II s-wave superconductor and 
then use magnetic flux tubes generated by an applied 
magnetic field to proliferate the vortices, which would 
then contain the Majorana states. In this work we seek to 
understand exactly when it is sufficient to only consider 
the surface physics, and when one must include the bulk 
physics. Very interesting work in this general direction 
is discussed in Ref. 15 where the role of chemical po- 
tential in the stability of the Majorana vortex modes is 
discussed for a topological insulator whose entire bulk 
has become superconducting. Here, instead, we focus 
only on the proximity effect scenario and the effects of 
the applied magnetic flux necessary to create a field of 
vortices. 

The manuscript is organized in the following manner: 
In Section |llj we detail the topological insulator Hamil- 
tonian utilized in this work. In Section|lIl] we review the 
physics resulting from the addition of very thin mag- 
netic flux lines in 3D topological insulators. In particu- 
lar, we review how a magnetic flux line connect the sur- 
faces of 3D topological insulators in which it enters and 
exits with a line of low-energy modes. In Section IVj we 



extend our analysis from the addition of magnetic flux 



lines in 3D topological insulators to topological insu- 
lators with s-wave superconducting pairing on the top 
and bottom surfaces. In this system, we discuss two dif- 
ferent physical regimes delineated by the spread of the 
magnetic flux as it penetrates the heterostructure. In the 
first physical regime, we study the behavior of the topo- 
logical insulator - superconductor heterostructure when 
the spread of the magnetic flux lines inserted into the 
system are limited in spatial extent to a size on the order 
of the lattice constant. This leads to the removal of the 
zero energy Majorana state from the system as the sur- 
face bound states may now tunnel along the magnetic 
flux tube and annihilate the states on the other surface. 
In the second regime, we study the case when the spread 
of the magnetic flux line has a much wider spatial ex- 
tent. In this situation, the Majorana fermions become 
localized at the interface between the topological insula- 
tor and proximity coupled superconductor and the bulk 
remains inert so that only the surface physics need be 
considered. 
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FIG. 1. (a) Schematic of a cylindrical 3D topological insulator 
with a hole drilled through the center. The blue line repre- 
sents a flux tube threaded through the cylindrical hole, (b) 
Schematic of a heterostructure of a topological insulator thin- 
film sandwhiched between two s-wave superconductors. The 
thin blue lines represent h/2e flux tubes which generate vor- 
tices in the superconductor layers. 



II. MODEL HAMILTONIAN 

In order to capture the essential physics of the prob- 
lem, we use a minimal bulk model for a 3D topological 
insulator which consists of a gapped Dirac Hamiltonian 

p 

= £4 (d a (p)T a + M(p)T°)c p . (1) 
p 

where a = 1, 2, 3, F 1 = t x ® a a , T° = t z <Ei I, a a is spin, 
r a is an orbital degree of freedom representing orbitals 
A, B, and c p = (c p Af c?Ai c pBt c pB i) T ■ In this work, to 
illustrate the salient physics, we will use both a contin- 
uum description with 

d a (p) = hv FPa , M(p) = m - {l/2)bp 2 (2) 

and a lattice description with 

d a (p) = (hv F /a)sin(p a a), (3) 

and 

M(p) =(b/a 2 ) (cos(p x a) + cos(p y a) + cos(p z a)) 
— 3b/ a 2 + m 

where VF,rn,b are material parameters and a is the lat- 
tice constant. These material parameters may be ad- 
justed using the previously tabulated constants based 
on DFT calculation !] 16 ! 17 ! to fit many of the most common 
3D topological insulators. Here, to simplify the nota- 
tion, we will set a and hv p equal to unity in the remain- 
der of the work unless otherwise noted. This model has 
time-reversal symmetry with T — I <g> ia y K where K is 
complex conjugation. For generic values of m ^ the 



system is a gapped insulator and we focus on the low- 
energy regime when m ~ 0.. Assuming translation sym- 
metry, the energy spectrum of the continuum model is 
E± = ± \Jp 2 + (m — (l/2)bp 2 ) 2 with each band doubly 
degenerate. As a convention, which is consistent with 
the behavior in canonical topological insulators such as 
Bi2Se3, we choose b > and as a result the trivial (topo- 
logical) insulator state occurs when m/b < (m/b > 0). 
In the following, when we refer to a topological insu- 
lator state, we are referring to a state described by the 
model Eq. [T|with m > and b > 0. 



III. MAGNETIC FLUX LINES IN 3D TOPOLOGICAL 
INSULATORS 

The physics of thin flux lines in the bulk of a topo- 
logical insulator was originally considered in Refs. [131 
and HU Let us begin with an infinite solid cylinder of 
3D topological insulator whose length is placed along 
the z-direction with a cylindrical hole drilled through 
the center as seen in Fig. [l|a). We take the inner and 
outer radii to be Rj, Ro respectively. Due to the char- 
acteristic property of time-reversal invariant topological 
insulators, there are low-energy modes bound to the in- 
ner and outer cylindrical surfaces. However, the surface 
fermions have a 7r-Berry phase when a particle winds 
around the Fermi surface. This leads to a condition that 
there will be no exact zero modes in the surface energy 
spectrum on the inner or outer surfaces, so long as we 
consider a cylinder of finite radius. To be clear, sur- 
face electrons that travel around the azimuthal direction 
on the inner or outer surfaces pick up a tt -Berry phase 
leading to effective anti-periodic boundary conditions 
which shifts the zero-momentum Fourier mode away 
from zero energy. To recover the exact zero-modes we 
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must twist the boundary conditions back to being peri- 
odic. This is accomplished by threading 7r-flux {4>q/2 = 
h/2e) through the hole drilled in the cylinder SEI 

In order to be concrete about the behavior of these 
zero energy modes, we begin with the continuum model 
Hamiltonian for a topological insulator introduced in 
Eq. FT] and assume the cylindrical with flux <f> in the unit 
of nfe threaded through the interior hole. Keeping only 
the linear terms in p we get 



H =niT z ® I + [p x - eA x )r : 
+ (Py 



eA y )r x <g> o y +p z r x 



(5) 



In Eq. (jsjl, the vector potential is A = ^ %{—y% + xy). As 
we have now added the magnetic flux into our Hamilto- 
nian, momenta p x and p y are not no longer good quan- 
tum numbers. It is important to note that The Hamil- 
tonian should be solved in real space where the mo- 
mentum operators are represented as p x ^ — id/dx and 
p y ^r — id/dy. First we will consider the case p z = and 
solve for the zero-mode eigenstates. Converting from 
Cartesian coordinates (x, y) to polar coordinates (r, 9) 
the Hamiltonian then becomes 



H, 



linear 



Po 




(6) 



As we are searching for the zero energy modes in the 
system, we must solve the eigenvalue problem for this 
matrix. For the energies we obtain 



Ri 



(7) 



with corresponding eigenstates 

e _ /« 7 m ( r ') dr ' 



4>i 



(0,-ze l ( £+1 ) 9 ,e l£e ,0) T , (8) 



ay2r 

where a is a normalization coefficient defined as 

poo 

2 n / -2 m(r )dr , 

a = 2ir / e Ja i dr. 



(9) 



In Eq. [7] we note that t, an integer which represents an 
angular momentum quantum number though it should 
be noted that for the eigenstates, different components 
may possess different angular momenta. Thus, if the 
flux is (j> = 1/2 + n for all integers n, then the resultant 
wavefunctions, ^0+ ar *d i/jq- f° r t — n, have zero energy. 



Following the procedure outlined in Ref . we now 
imagine adiabatically shrinking the radius Rj — >■ a and 
consider a single lattice plaquette as the hole drilled 
through the center of the cylindrical topological insu- 
lator. Thus, 7r-flux threaded through a line of single- 
plaquettes produces zero modes localized on the line of 
plaquettes along the length of the cylinder on the in- 
ner and the outer boundary. If we turn on p z we will 
find a Kramers' pair of propagating modes on the inner 
and outer surfaces with disperse linearly in p z . As we 
have time reversal invariance, we expect modes prop- 
agating in both directions with opposite spin polariza- 
tions. Therefore, a 7r-flux line confined to a hole, even in 
the limit where the hole is reduced to the size of a single 
plaquette, in a topological insulator will trap a single- 
pair of gapless counter propagating modes akin to the 
ID holographic edge state found in a 2D quantum spin 
Hall system. This is referred to as the wormhole effecP 
and it demonstrates that the bulk of a topological insula- 
tor is not generically inert when in the presence of mag- 
netic flux. 



IV. 3D TOPOLOGICAL INSULATOR - 
SUPERCONDUCTOR HETEROSTRUCTURE 

Having shown that the presence of magnetic flux in a 
3D topological insulator forces one to consider the pres- 
ence of a non-inert bulk, we proceed to understanding 
the effects of coupling type-II s-wave superconductors 
to the surfaces of a 3D topological insulator. We will 
want to add a sufficient amount of magnetic flux to gen- 
erate vortices yet not so much as to necessitate the con- 
sideration of the interactions or quasi-particle tunneling 
between vortices. With proximity coupling to an s-wave 
superconductor, and in the presence of a magnetic field 
B, we must use the Bogoliubov-de-Gennes (BdG) mean- 
field description of our Hamiltonian: 



H 



BdG 



H D (p-eA) 
At 



A 

-HU-V 



eA) 



(10) 



with V x A = B, A = A (x)I ® ia y , and * p = 



p l. pi . We consider a heterostructure geometry 

with a thin-film of topological insulator sandwiched 
along the z-direction between two s-wave superconduc- 
tors as shown in Fig. [ljb). We model the physics of the 
superconductors by inserting an induced s-wave pair- 
ing term into the BdG Hamiltonian which penetrates 
into the topological insulator film. Without the pres- 
ence of a magnetic field, this implies that Ao (x) = Aq(z). 
That is, we assume that that the superconducting prox- 
imity pairing is homogenous in the xy-plane if no vor- 
tices are present. 

With this in mind, we can proceed with our analysis of 
the effects of the pairing. It is important to note that the 



4 



pairing, if weaker than the energy scale of the bulk insu- 
lating gap, will not affect the gapped bulk states of the 
topological insulator. However, it will affect the metallic 
surface states which are susceptible to a superconduct- 
ing pairing potential. The effective surface BdG Hamil- 
tonian is 

rr(surf) _l^ $ t/ Px<J V ~ PyCT X A ia y \ , 

BdG - 2 L*p^ -A*i<jy -p x uy - Py a x ) p 

(ii) 

where $ p = [c^ c p ^ cl p ^ c'-pij and A represents 

the effective pairing potential felt by the surface states. 
This Hamiltonian has a gapped energy spectrum E± = 
± yjp 2 + |A | 2 . Thus, a non-zero proximity coupling in- 
duces a gap in the topological surface states. Previous 
work has shown that a vortex induced on the proximity- 
coupled surface traps a Majorana bound state. This is 
shown by solving Eq. [TTlwith a vortex present, which is 
inserted by winding the superconducting order param- 
eter A as 

A - A (r)e w ^ (12) 

where 6(r) is the polar angled 

While it may be possible to find Majorana states at the 
center of vortices in 3D topological insulator - s-wave 
superconductor heterostructures, there is not a standard 
prescription of how to generate such vortices. We con- 
sider the simplest possible route and apply a uniform 
external magnetic field perpendicular to the heterostruc- 
ture. Physically, we must apply a large enough mag- 
netic field to generate vortices, but small enough that 
the vortex density is low as mentioned above. As the 
topological insulator film is not inert to the addition of 
flux, we must be careful to account for the effects of 
the applied magnetic flux to ensure it does not spoil 
the bound state structure. For simplicity, we consider 
only the one (analytic results) and two-vortex (numeri- 
cal results) problems noting that the one and two vortex 
problems are essentially equivalent, and only differ be- 
cause of the choice of boundary conditions. If our sys- 
tem contains periodic boundary conditions in the x and 
y directions then we must have an even number of vor- 
tices; this is the situation we consider in our numerics. 
If we choose open boundary conditions in x and y, then 
a single-vortex in the bulk implies the existence of an- 
other vortex at the boundary or at infinity; this is the 
case for our analytic results. We assume that a magnetic 
flux of only one <fio quantum, parallel to the z -direction, 
penetrates the superconductors. In our heterostructure 
the superconductors on top and bottom would, in prin- 
ciple, dynamically generate two vortices in each layer. 
We will assume that the induced vortices (that we put 
in by hand) are well-separated enough so that they do 
not influence each other and that the positions of the 
vortices on the top and bottom surfaces share the same 
(x, y) position for simplicity. 



Inside the superconductor the penetrating magnetic 
field satisfies the London equatiorPO 

B(r) - A 2 V 2 B(r) = ^(r) (13) 

near a vortex positioned at the origin with penetration 
depth A. The solution for B(r) in the superconductor is 

B(t) = z^K (r/X) (14) 

where K n (x) are modified Bessel functions of the 
second-kind. The flux within the disk of the radius r 
is 

0(r) = (l/2-(r/2A)ifi(r/A)). (15) 

We want to model the effects of B(r) in the entire het- 
erostructure including the topological insulator but this 
not easy to account for. We instead opt for a more 
phenomenological approach to capture the qualitative 
physics. Once the flux leaves the superconducting 
layers and enters the topological insulator film it will 
spread out. For our purposes, we will consider a model 
where A varies with the depth in the heterostructure as 
A = \(z) and study how the vortex physics changes 
with X(z). If the film is thin, the flux will not have suf- 
ficient distance to spread before it must re-enter the top 
superconducting layer and thus modeling the insulator 
layer as having a finite penetration depth (larger than 
that of the superconductor) is not unreasonable. In or- 
der to understand the appropriate physics, it is natural, 
in the context of this problem to consider two separate 
limits associated with the amount of magnetic flux pen- 
etration into the topological insulator. In the first limit, 
we wish to examine the "thin-flux" limit in which the 
flux which penetrates the topological insulator does not 
spread out very far in the topological insulator before 
reentering the other superconducting layer. In the sec- 
ond limit, we examine the case in which the magnetic 
flux spreads out widely in the topological insulator be- 
fore it mush re-enter the other superconducting layer. 
These two limits can be considered analytically while 
we provide numerical calculations which capture the in- 
terpolation between these cases. 

A. Thin-Flux Limit (A ~ a) 

We begin from the limit of two well-separated, thin 
flux tubes of flux <fio/2 where the flux tubes are each 
confined within single plaquettes i.e. A <C a. When the 
proximity pairing potential vanishes, the system will ex- 
hibit gapless modes propagating on each of the thin flux 
tubes (ignoring finite size splitting due to hybridization 
with the second vortex). The resulting gapless theory 
of a single tube is simple to understand using the re- 
sults from the previous section. There we solved Eq. [T] 
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at p z — with a 7r-flux tube through a single plaque- 
tte to obtain the two zero-mode solutions V'o+iV'o- ( a 
Kramers' pair) localized on the flux. Then we can use 
k • P perturbation theory and treat p z as a perturbation to 
obtain, in the basis of tpo+^o-r the low-energy Hamil- 
tonian Hfi ux -u ne = p z a x . This Hamiltonian is iden- 
tical to the edge Hamiltonian of a quantum spin Hall 
edge state, as mentioned earlier. If we begin to increase 
A which corresponds to allowing the magnetic flux to 
spread uniformly in the z-direction, i.e. we move away 
from the wormhole limit, this applies a perturbation to 
the gapless flux-line Hamiltonian. Using perturbation 
theory, we find that 



of freedom are localized near the flux line so we can use 



Hfiux-une = Pz<J x + m x (X)a v + m y (X)a z 



(16) 



where the mass term, m,-, is monotonically increasing 
as A increases. This Hamiltonian has a gapped energy 

spectrum E± = ± ^Jp 2 z + m 2 + rn% which is expected 

since time-reversal is broken and the flux-line Kramers' 
degeneracy at p z = is lifted. Note that we are not in- 
creasing the amount of flux, only the region over which 
it spreads. If the 7r-flux tube is larger than one plaque- 
tte then some bonds in the lattice model will necessar- 
ily have phase factors which have imaginary contribu- 
tions to the Hamiltonian, regardless of the gauge choice, 
which break the time-reversal symmetry of the system. 
Our perturbation theory analysis is approximately valid 



until the induced gap E m = y m 2 + m 2 approaches the 

bulk mass gap m. To estimate the size of the induced 
gap, Em, caused by the spreading of the flux in the topo- 
logical insulator (<fi(r)) in Eq. 15 we apply first-order 
perturbation theory: 

E M = (^ 0+ | AH | 7/W) = -<</>o- I AiJ | Vo-), (17) 



where AH = H Hnear (<j>(r)) - H u 



1/2). Using 



the previously obtained expression for Hu near in Eq. pi 
the first order approximation for Em in the continuum 
limit is 



E 



7T 

Xa 2 



K 1 (r/X)e 



dr 



(18) 



We refer to this flux regime as the 'thin-flux limit' i.e. 
the regime in which we can consider the low-energy 
states as those originating from the gapped wormhole 
modes. With the effects of the magnetic flux accounted 
for, we now turn on the superconductor proximity effect 
in the thin-flux limit. There will be an induced super- 
conducting pairing potential that is z-dependent and, 
when flux and the corresponding vortices are present, 
the superconducting pairing takes on an x and y depen- 
dence as well. Before we get to the situation where Ao 
is only non-vanishing near the top and bottom surfaces, 
let us consider an induced Ao which is homogenous in 
the z-direction over the entire topological insulator. In 
the thin-flux limit the only low-energy metallic degrees 



our effective flux-line Hamiltonian from Eq. 16 to form 
a BdG Hamiltonian for the low-energy degrees of free- 
dom: 



H 



(BdG) 
flux— line 



(P) 



1 / H 



flux— line 



(p) 



iA a y 



-iA* aV -H* flux _ Une (-p) 



which has an energy spectrum with four non- 
degenerate bands 



±E+ 



|A |±£ 



M 



(19) 



This spectrum is gapped unless | A 1 = \E M \ - Now let us 
consider more realistic conditions where the thin-film is 
too thick to become entirely superconducting, and the 
proximity induced pairing depends on z. Specifically, 
the superconducting pairing potential decays as we 
move away from the surfaces towards the interior the 
topological insulator film. In this case the system has a 
time-reversal symmetry breaking mass Em which is ho- 
mogenous in the z-direction (as per our phenomenolog- 
ical model) and a superconducting mass, | Ao(z)|, which 
is z-dependent. From standard ID Dirac physics^^ 
this model will exhibit localized, zero-energy Majorana 
bound states on mass domain walls in the z-direction 
along the flux line i.e. the places where |Ao(z)| = |-Em|- 
As the thickness of the flux increases so does Em and 
the domain-walls along the vortex line get pushed to- 
ward the surface. 

This perturbation theory is valid as long as Em <C m. 
If we want to be able to carry out a full interpolation be- 
tween the thin-flux limit and the thick flux limit (to be 
discussed in the next section) we must rely on a numer- 
ical calculation. We show the results of such a calcula- 
tion in Fig. [2] We used a Lanczos exact-diagonalization 
algorithm to solve for the zero-modes of a full 3D lat- 
tice model. The vortices and proximity effect associated 
with the superconducting regions were non-dynamical 
and included in the mean-field limit following, for ex- 
ample, Ref . 20 We present the details of our numerical 
calculations in Appendix A. As illustrated in Fig. |2k-f, 
when we allow the flux to spread in the topological in- 
sulator film, i.e. as A increases, the domain-wall bound 
states, which begin in the interior of the topological in- 
sulator, move outward toward the surface. As discussed 
in the previous paragraph this can be understood by 
noting that as A increases Em increases and the position 
of the mass-domain wall moves toward the surface. As 
the flux becomes thicker the bound states become more 
localized on the surfaces at the points around which the 
superconducting order parameter winds due to the flux. 
Since there are two domain-walls on each flux tube the 
pair of bound states will hybridize and lie higher than 
zero energy as shown in Fig. [2jg. As A increases the hy- 
bridization decreases which rapidly drives the states to- 
wards zero energy. 
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FIG. 2. Probability distribution for the lowest energy states corresponding to different superconducting penetration depths, A 
with: (a) A = 0.001, (b) A = 0.2, (c) A = 0.3, (d) A = 0.4, (e) A = 0.5, (f) A = 1.5. Note that A is in units of the lattice constant 
a and the entire flux is spread out in a region with a radius of roughly 5A. The pairing potential A(z) decays from 0.5 On the 
surface to within 5 layers. As A is increased we see that the states move from being delocalized along the flux tube penetrating 
the bulk of the topological insulator to being pinned at the surface. The inset shows a schematic of the spatial variation of the 
superconducting mass and the time-reversal breaking mass associated with the magnetization as A varies. (g)The energies of the 
lowest energy states as a function of A. There are clear zero modes forming as A — > oo. Our numerical results show that in A < a 
regime, the lowest energy linearly decays as the height of the sample increases; in A 2> a regime, the lowest energy exponentially 
decays as the height of the sample increases. 



B. Thick-Flux Limit (A > a) 



From Fig|2j we see that in the extreme thin-flux limit, 
the Majorana modes will penetrate into the bulk and hy- 
bridize with the states on the other surface and annihi- 
late. Fortunately, as indicated by our analytic pertur- 
bation theory, and numeric lattice model calculations, 
A does not have to be very large before we move from 
the wormhole effect /thin-flux limit so that the vortex 
modes are tightly bound to the surface at zero energy. 
The wormhole effect simply generates a region in the 
bulk with a mini-gap across which the Majorana states 
can tunnel to the opposite surface. While the bulk of a 
topological insulator is not inert to flux insertion, as long 
as the flux is not tightly bound to a region on the order 
of a lattice plaquette, the Majorana states will have dif- 
ficulty tunneling between the top and bottom surfaces 
and will be well-localized in the surface vortex cores. 

Once the flux is thick enough to restore the bulk gap 
entirely, we can consider the explicit Majorana bound 
state solution in the presence of the magnetic flux from 
the surface Hamiltonian alone. We consider the BdG 



The surface Dirac Hamiltonian is 



surface-state Hamiltonian in Eq. 11 with non-zero vor- 
tex winding and magnetic flux. We focus on the neigh- 
borhood of a single vortex and solve the problem for 
generic flux and order parameter profiles in the contin- 
uum limit. We begin by assuming we have a vortex at 
the origin generated by a magnetic flux given by Eq. 14 



H 



H(p, 



ia y Aoe 



-ia v A*e l " -H*(-p, 



(20) 



-if.) 



hvj. 



r 



-id g +4> 



-ida+4> 



where we have changed to polar coordinates and have 
implemented a non-zero vector potential Ag — H<f>(r) / er 
where 4>(r) is given in Eq. 15 This Hamiltonian has two 
eigenstates with zero energy: 



1 f rf A(r') <t>(r') \ Hr , 

= -e~ Jo {-^ +—) dr 



|V> 2 ) = -e~ J ° 1^+ — )* 



/ e~ i9 \ 



V o J 



(21) 



(22) 



where /3 is a normalization coefficient. If we turn off the 
magnetic field i.e. for cf> = we are in the Fu-Kane limit 
with very diffuse flux and then only the state \ipi) is nor- 
malizable which matches their result 9 . For a generic flux 
profile there will still only be one zero-mode solution 
that satisfies the boundary conditions and it will be a 
linear combination | t/> ) = a i(0)|'0i) + ^(fi)]^) ■ The co- 
efficients ai (<f) , ci2 ((/)) control the spin composition of the 
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zero-mode and depend on not only the detailed bound- 
ary conditions but also on the short-distance physics of 
the vortex structure. 

Since the coefficients a\ (</>), ai (<fi) depend on the de- 
tails of the system we numerically calculate them. By 
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solving the eigenvalue problem of H^Iq from Eq. 
a zero mode Majorana bound state exists in the core o 
a vortex. The ratio of spin up to spin down (|ei2|/|ai|) 
for a single zero mode is shown in Fig. [3] This ratio 
describes the mixing between the two allowed zero en- 
ergy modes when finite flux is present. We find that at 
A ~ 0.1, where the magnetic flux starts to spread over 
more that one plaquette, the ratio of spin up and down 
starts to decrease rapidly. As A — > oo, 1 02 1 / 1 o-x | — > 0. In 
this limit, this is a Fu-Kane Majorana bound state 2 -^ pos- 
sessing a single species of spin. As A — > 0, |a2|/|ai| — » 1. 
This limit is the thin-flux limit, for which we see that 
the zero modes have the same portions of spin up and 
down. This is due to the fact that in Eq.[l6]the Hamilto- 
nian Hfi U x-Une contains equal portions of spin up and 
down, and, therefore, the zero modes in this limit have 
the same portions of spin up and down. In short, as A 
decreases we move from the thick-flux to thin-flux limits 
and |oi(^)|(|o2(^)|) monotonically decreases (increases). 
Thus, we find that in the limit where the flux does not 
affect the bulk physics the effective magnetic field only 
acts to change the spin composition of the zero-energy 
vortex core state. 




log 10 X 



FIG. 3. Ratio of spin up to spin down composition of a Majo- 
rana bound state versus flux penetration depth A. The differ- 
ent traces represent the use of different superconducting pair- 
ing potential strengths (A). In the thick-flux limit, a larger A 
corresponds smaller ratio of spin up and down leading to spin 
polarized Majorana bound states. 



Beyond understanding the spatial extent of the flux 
spread on the Majorana states, we wish to look at chang- 
ing A, which changes the extent of the Majorana bound 
state. We fix the penetration depth of the flux to lie in 
the thick-flux regime and consider the probability dis- 
tribution of the zero modes. In the thick flux limit, we 
find that the spin down dominates the composition of 
the Majorana bound states. Therefore, zero-mode wave- 
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FIG. 4. Probability distribution of a Majorana bound state in 
the in the topological insulator / s-wave superconductor het- 
erostructure where we have used realistic materials parameter. 
For the topological insulator we have used the parameters of 
Bi2Se3 and for the superconducting films on the top and bot- 
tom surfaces we have used the material parameters of a nio- 
bium superconductor with a single vortex. 



function, described in Eq. 21 has a probability distribu- 
tion 



P(r) = (V>1 | fa) = JL e -*Ar/Hv F -K (r/A)_ (23) 



The decay length of the probability distribution approxi- 
mately equals Hvf/A for A 3> Hvf/ A. The reason is that 

as x -> 0, K Q (x) In a; so then P(r) ~ e -2Ar/hv P 

which is the Fu-Kane result 9 . However, for A < hvp/ A, 
we have to consider the probability distribution in Eq. 
[23]directry to find the width of the zero modes. 

It is important to solidify these results by discussing 
our results within a realistic context. Therefore, we esti- 
mate the width of Majorana bound states with real phys- 
ical parameters. According to Ref .|l7l for the Bi2Se3 fam- 
ily of topological insulators Hvf ~ ^eVA. For the super- 
conducting top and bottom layers we use the pairing 
potential of the type-II superconductor niobium which 
approximately equals 1 meV. We assume that the prox- 
imity effect, which results in a pairing potential on the 
surface of the topological insulator has an induced pair- 
ing also of the order of 1 meV which is the best-case sce- 
nario. Using these parameters, the quantity hvp/A is 
about 400 nm, which is much larger than the London 
penetration depth of 40 nm 22 . This combination of pa- 
rameters allows us to now plot the probability distribu- 
tion from Eq. 23 directly as shown in Fig. E] We find 
that in this case, the width of a Majorana bound state 
is of the same order of magnitude as that of the pen- 
etration depth. In general, the Majorana bound states 
will have the same characteristics as shown in Fig. [4] for 
A < hvp/A. In most experimental regimes the physics 
will be deep in the thick-flux region where only the sur- 
face physics is important. One notable exception would 
be experiments where vortex pinning sites are artifi- 
cially created (e.g. by drilling through the heterostruc- 
ture). In this case there is a possibility that the holes 
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inside the topological insulator could trap a sizable frac- 
tion of a 7r-flux quantum which will lead to a mini-gap 
region across which the majorana bound states can tun- 
nel. 



V. CONCLUSION 

In conclusion, we have illustrated the effects of mag- 
netic flux in topological insulator / superconductor het- 
erostructures in two different regimes. In the thin-flux 
limit the Majorana fermion bound states can be destabi- 
lized through hybridization with low-energy bulk states 
localized near the thin flux line. These effects will 
be more pronounced in thin topological insulator films 
with minimal flux spreading, or in samples where vor- 
tex pinning sites are produced by drilling holes through 
the heterostructures. If such holes continue to trap ap- 
proximately h/2e flux throughout the topological insu- 
lator film then effects of flux in the bulk must be care- 
fully considered. The opposite regime, where the vortex 
core does not feel much effective flux is likely the phys- 
ical regime of most experiments. In this case we can ig- 
nore the bulk effects and focus only on the surface. The 
flux in this regime simply acts to change the spin con- 
tent of the vortex zero mode and does not affect its en- 
ergy or stability. Furthermore, when we consider the pa- 
rameters corresponding to real topological insulator / s- 
wave superconductor heterostructure with Bi2Se3 as the 
topological insulator and with niobium superconductor 
layers, we find that a Majorana fermion trapped in the 
magnetic flux is stable with a spatial extent of around 40 



nm. 
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Appendix A: Numerical Calculations for the Bulk 
Hamiltonian 

We construct the Hamiltonian of a strong topologi- 
cal insulator sandwiched between two s-wave super- 
conductors with vortices. The lattice Dirac model we 
use for the 3D topological insulator is 

H D = J2 H(m)4c r + £ H(S)4c r+s (Al) 

r r,<5 

hi \ ql tt/x\ bT + iA6-f 
H (m) — m — So, H (o) = , 

where r indicates the position of the lattice, 5(= 
±ax, ±ay, ±az) indicates the nearest neighbor hopping, 
and r = Fix + T^y + ^3%. Consider the interface be- 
tween a strong topological insulator and an s-wave type 
II superconductor with an even number of vortices. The 
proximity effect leads to the pairing potential of the 
superconductor leaking into the topological insulator. 
Near the interface of the topological insulator we can 
write down the 8x8 BCS-type lattice BdG Hamiltonian 



/ t W H(m) A I <g> ia y e^ \ ( c r 

H BdG ~ ^ ( 4 c r ) [ ^ ^ *( p) _ H . (m) ) { c t 



E 

r,5 



+ ^ (4Cr) ( -2r(*)e*f/T»A«.* )U + J' (A2) 



where A is the vector potential is coming from the magnetic field B described by the conventional London equation 



with penetration depth A, as in Eq. 14 with the cores of the vortices at rj. The phase </>(r) acts as an additional 
"gauge field" coupled to the quasiparticles. With vortices, <j>(r) is not a pure gauge: V x V0(r) = 277^^^ 5(r — 
rj). For the numerical calculation, we want to cancel out the phase <fi(r) in the order parameter to speed-up the 
simulation. Although 4>{r) is not a pure gauge, by performing a "bipartite" singular gauge transformation 2 ^ the 
phase is successfully moved to the diagonal terms of the Hamiltonian. Here, in our main numerical calculation, we 
consider two vortices located at r A and r B respectively. A "bipartite" singular gauge transformation is c r = c' r e l ^ A ^ r ^ 
for the particle part and c r = cJ,e^ B ^ r ^ for the hole part, where V x V4>a{y) — 2nzS(r — r A ) and V x V0s(r) = 
2nzd(r — r B ). This gauge transformation avoids a multi-valued problem so that the integral f* +d \7<fi A / B (r) ■ dl = 
( l ) A/B{ r + 8) ~~ 4>a/b{ y ) is path-independent up to 2itn, which does not affect the probability distributions of the 
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eigenstates from the numerical calculation. The Hamiltonian we use in the simulation becomes 

HbcIG = ^2 1 



C r Cr 



where 



H(m) A I <S> i<J y \ ( c r 
-A I ®io y -H*(m) J I cj 



r,<5 



^ + \v<j> A/B {v)~ e -A{v)]-d\ 



1 r +,5 r r (;/-?/ A/i3 h (a;-^/ 5 ) (y~y B/A ) , (x~x B / A ) 



The size of the topological insulator is (n x — 1) X (n y — 
1) x (n z — 1). In the numerical calculation, we typically 
used n x — 28, n y — 20, and n z = 24 with the lattice 
constant a = 1, also with open boundary conditions in 
all directions. We are modeling a thin-film of topologi- 
cal insulator sandwiched along the z-direction between 
two s-wave superconductors. Because of the proximity 
effect, we assume A = A (z) smoothly decays away 
from the top and the bottom surfaces and vanishes in 
the middle region. In the xy-plane, let the center be the 
origin (0,0). The positions of the two vortices are set at 
(ncc/4, 0) and (— nx/A, 0). Also in the thin film limit, we 
assume phenomenologically that the penetration depth 
A is independent of z. We choose m = 1.5, 6=1, and 
A = Hvp = 1. This is a strong topological insulator 
phase with topological invariants (1; 111). Finally, after 
solving the eigenvalue problem to find the lowest en- 
ergy modes of the Hamiltonian in Eq. A3 the probabil- 
ity distributions of the lowest energy modes with vary- 
ing A are shown in Fig. [2] 



Hamiltonian can be written explicitly 

HsdG = 

i(M (p)cr 2 ® T° + smp x l (g) T 1 + sinp y (T z ® T 2 

+ sinp z I ® T 3 + A R a x <g> I <g> ia y — A/<r y <g> I <g> ia y ), 

(Bl) 

where An and Aj are the real and imaginary parts of the 
order parameter. We start by finding the zero modes on 
the surface of a strong topological insulator. Therefore, 
to find the surface/domain-wall zero modes we need 
to have H z \ijj) = (M(p)a z ® T° + sinpj <8> r 3 )|-0) = 0. 
Qualitatively we can assume near the surface for z > 
0, m is positive, and for z < 0, m is negative. Hence, 
p z is not a good quantum number. For the low energy 
physics, i.e. when we focus around m ~ 0, k ~ we can 
safely take the continuum limit so that sinp z — > —id/dz. 
The wavefunction of the surface state is proportional to 
e~ Jo m ( 2 ) dz . There are four zero modes solutions: 



Appendix B: Numerical Calculations for the Surface 
Hamiltonian 



In the thick-flux limit, the zero modes are pushed 
to the surface of the topological insulator. Therefore, 
the surface Hamiltonian can adequately describe the 
physics of the zero modes. Solving an eigenvalue prob- 
lem of the 2D Hamiltonian allows us to consider a larger 
size of the system. In the following, we will derive the 
surface Hamiltonian from the bulk BdG Hamiltonian in 



Eq. A3 with vanishing magnetic field, and transform 
it to position space for the simulation. The bulk BdG 



|pt) =F(z)(-l,0,i,0) T , 
| fct) =F(*)(-l,0,-i,0) T , 

p~ Jo m ( z ') dz ' 

F(z) =- 



N 



\pi) =F(z)(0,l,0,z) T 
| hi) =F(2)(0,1,0,-*) T , 

(B2) 



where N is a normalization coefficient, p and h indicate 
particle and hole parts respectively, and t, I are associ- 
ated with spin up and down. Thus, the projection of 
the bulk HsdG to these four modes is an effective sur- 



face Hamiltonian, which was written in Eq. 11 We note 
that if the boundary condition changes (M — > —M), 
the surface Hamiltonian is the same in the similar basis 



of (\E'^, ^i, i&l, ^T) while the basis wavef unctions of the 
four zero modes change. For the top and the bottom sur- 
faces, the physics can be described by the same surface 
Hamiltonian. For convenience we use a simple lattice 
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regularization for the numerical calculation of the 2D regularization we use a 2D lattice Dirac model tuned to 
^BdG m position space with two vortices. For the lattice the critical point: 



n surf _ \p / + \ ( h{m) Aia y \ ( c r 
BdG - 2^\ c r Cr ) I -Ma y ~h*{m) )\c\ 



E 



t , / ft( £ )e«/r r+ '[v^(r)- lfe A W ].«a o 

Cr Cr 1 [ o -/ l *(e) e - i / r r+e [ v ^ B «-^ A ( r )]" il 11 J 1 ' 




/i(m) = ?n(T 2 , /i(e) 



cr 2 + «e ■ 7 



where 7 = (<r„, — cr-,.) and the nearest neighbor hopping is described by e = ±i and ±y. We set m = —2 to have a 
gapless two-component, 2D Dirac cone Hamiltonian when A vanishes. To avoid boundary effects, we use periodic 
boundary conditions in (x, y). We used a size of the surface of n x x n y = 120 x 60. The positions of the two vortices 
are (±30, 0). The surface Hamiltonian calculation shows the ratio of spin up and down in Fig. [3] The use of the lattice 
Hamiltonian is only valid if we are looking for low-energy properties of the spectrum, and for example, it does not 
satisfy the same symmetry properties under time-reversal that a true surface state Hamilotnian of a 3D topological 
insulator would. 
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